ecospat优化可视化

### 另外的教程做可视化：

library(knitr)
library(spThin)
library(rgeos)
library(sp)
library(maptools)
library(raster)
library(ecospat)

occ.points <- read.table("points.txt", sep = "\t", header = TRUE)
occ.points.thin <- thin(occ.points, verbose = FALSE, 
                        lat.col = "Latitude",
                        long.col = "Longitude",
                        spec.col = "Scientific.name",
                        thin.par = 10,
                        reps = 1, 
                        write.files = FALSE,
                        write.log.file = FALSE,
                        locs.thinned.list.return = TRUE)

data(wrld_simpl)plot(wrld_simpl)points(occ.points.thin[[1]], col = "purple", pch = 20, cex = 0.7)

## 根据分布区设置不同的颜色：## 这一部分可以用arcgis替代，作图效果更好；
n.groups <- 5 
g.names <- c("Western North America",
             "Eastern North America",
             "Native",
             "Western Australia",
             "Eastern Australia")

g.codenames <- c("WestNorAme", "EastNorAme", "Native", "WestAus", "EastAus")
g.colors <- c("cyan", "darkblue", "red", "green", "darkgreen")

## 背景定义：
buffer.size <- 1
mcp <- function (xy) {
  xy <- as.data.frame(coordinates(xy))
  coords.t <- chull(xy[, 1], xy[, 2])
  xy.bord <- xy[coords.t, ]
  xy.bord <- rbind(xy.bord[nrow(xy.bord), ], xy.bord)
  return(SpatialPolygons(list(Polygons(list(Polygon(as.matrix(xy.bord))), 1))))
}

variables <- getData('worldclim', var='bio', res=10)
variables <- subset(variables, 1:19)

plot(variables)

# Union of the world map
lps <- getSpPPolygonsLabptSlots(wrld_simpl)
IDFourBins <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
world <- unionSpatialPolygons(wrld_simpl, IDFourBins)

# Empty objects
g.assign <- numeric(nrow(occ.points.thin[[1]]))
xy.mcp <- list()
back.env <- list()
spec.env <- list()
row.sp <- list()

# Plot map
plot(wrld_simpl)

# Loop
for(i in 1:n.groups) {

  # Define groups
  cut1 <- occ.points.thin[[1]][, 1] >= group.long[i]
  cut2 <- occ.points.thin[[1]][, 1] < group.long[i + 1]
  g.limit <- cut1 & cut2

  # Save row numbers per species
  row.sp[[i]] <- which(g.limit)
  g.assign[g.limit] <- g.names[i]

  # Background polygon
  mcp.occ <- mcp(as.matrix(occ.points.thin[[1]][g.limit, ]))
  xy.mcp.i <- gBuffer(mcp.occ, width = buffer.size)
  proj4string(xy.mcp.i) <- proj4string(world)
  xy.mcp[[i]] <- gIntersection(xy.mcp.i, world, byid=TRUE, drop_lower_td=TRUE)
  # Background environment
  back.env[[i]] <- na.exclude(do.call(rbind.data.frame, extract(variables, xy.mcp[[i]])))
  # Species environment
  spec.env[[i]] <- na.exclude(extract(variables, occ.points.thin[[1]][g.limit, ]))

  # Plot
  points(occ.points.thin[[1]][g.limit, ], col = g.colors[i],
         pch = 20, cex = 0.7)
  plot(xy.mcp[[i]], add = TRUE, border = g.colors[i], lwd = 2)

}


## 最终建模所使用的环境变量和数据：
# Occurrence points per group
g.occ.points <- cbind("Groups" = g.assign, occ.points.thin[[1]])
# Environmental values for the background 
all.back.env <- do.call(rbind.data.frame, back.env)
# Environmental values for the species occurrence points 
all.spec.env <- do.call(rbind.data.frame, spec.env)
# Environmental values all together
data.env <- rbind(all.spec.env, all.back.env) 
table(g.occ.points[, 1])

######################################################
################## niche comprision ################
####################################################
w <- c(rep(0, nrow(all.spec.env)), rep(1, nrow(all.back.env)))
# PCA of all environment
pca.cal <- dudi.pca(data.env, row.w = w, center = TRUE, 
                    scale = TRUE, scannf = FALSE, nf = 2)


# Rows in data corresponding to sp1
adtion <- cumsum(c(0, sapply(back.env, nrow)))
begnd <- nrow(all.spec.env)
# Empty list to save the results
scores.back <- list()
scores.spec <- list()

# Assigning the values 
for (i in 1:n.groups) {
  scores.spec[[i]] <- pca.cal$li[row.sp[[i]], ]
  pos <- (begnd[1] + adtion[i] + 1) : (begnd[1] + adtion[i + 1])
  scores.back[[i]] <- pca.cal$li[pos, ]  
}

total.scores.back <- do.call(rbind.data.frame, scores.back)

# Rows in data corresponding to sp1
adtion <- cumsum(c(0, sapply(back.env, nrow)))
begnd <- nrow(all.spec.env)
# Empty list to save the results

scores.back <- list()
scores.spec <- list()
# Assigning the values 
for (i in 1:n.groups) {
  scores.spec[[i]] <- pca.cal$li[row.sp[[i]], ]
  pos <- (begnd[1] + adtion[i] + 1) : (begnd[1] + adtion[i + 1])
  scores.back[[i]] <- pca.cal$li[pos, ] 
}


total.scores.back <- do.call(rbind.data.frame, scores.back)


R <- 100

z <- list()
for (i in 1:n.groups) {
  z[[i]] <- ecospat.grid.clim.dyn(total.scores.back,
                                  scores.back[[i]],
                                  scores.spec[[i]],
                                  R = R)
}

## 
# Once the number of interactions is de€ned, we can generate the values. Additionally,
# we calculate the partition of the non-overlapped niche, among niche un€lling,
# expansion and stability (see methods in the manuscript)
D <- matrix(nrow = n.groups, ncol = n.groups)
rownames(D) <- colnames(D) <- g.codenames
unfilling <- stability <- expansion <- sim <- D
for (i in 2:n.groups) {

  for (j in 1:(i - 1)) {

    x1 <- z[[i]]
    x2 <- z[[j]]

    # Niche overlap
    D[i, j] <- ecospat.niche.overlap (x1, x2, cor = TRUE)$D

    # Niche similarity 
    sim[i, j] <- ecospat.niche.similarity.test (x1, x2, rep,
                                                alternative = "greater")$p.D
    sim[j, i] <- ecospat.niche.similarity.test (x2, x1, rep,
                                                alternative = "greater")$p.D

    # Niche Expansion, Stability, and Unfilling
    index1 <- ecospat.niche.dyn.index (x1, x2, 
                                       intersection = NA)$dynamic.index.w
    index2 <- ecospat.niche.dyn.index (x2, x1,
                                       intersection = NA)$dynamic.index.w
    expansion[i, j] <- index1[1]
    stability[i, j] <- index1[2]
    unfilling[i, j] <- index1[3]
    expansion[j, i] <- index2[1]
    stability[j, i] <- index2[2]
    unfilling[j, i] <- index2[3]
  }
}

# D value:
kable(D, digits = 3, format = "markdown")


# Niche similarity null model (p-values):
kable(sim, digits = 3, format = "markdown")

## Niche Un€lling:
kable(unfillin diits = 3 format = "markdown")

## Niche Expansion:
kable(expansion, digits = 3, format = "markdown")

## Niche Stability:
kable(stability, digits = 3, format = "markdown")